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Abstract 

We consider polynomial systems of Prony type, appearing in many areas of 
mathematics. Their robust numerical solution is considered to be difficult, es¬ 
pecially in “near-colliding” situations. We consider a case when the structure of 
the system is a-priori fixed. We transform the nonlinear part of the Prony sys¬ 
tem into a Hankel-type polynomial system. Combining this representation with 
a recently discovered “decimation” technique, we present an algorithm which 
applies homotopy continuation to an appropriately chosen Hankel-type system 
as above. In this way, we are able to solve for the nonlinear variables of the 
original system with high accuracy when the data is perturbed. 


1. Introduction 

1.1. The Prony problem 

Consider the following approximate algebraic problem. 

Problem 1. Given (mo,... ,rhN-i) G C^, e ^ 0, s S N, a multiplicity vector 
D = {di,... ,ds) G N® with d := 2d y N, find complex numbers 

{zj, satisfying ^ 0 and \zj\ = 1 with {zj} pairwise dis¬ 

tinct, such that for some perturbation vector (cq, •. •, ew-i) S C'^ with \ek\ < £ 
we have 

s — 1 

mk +ek, k = 0,...,N -1. (1) 

i=i e=o 

^ ^ > 

:=mk 


Email address: batenkov@cs.teclmioii.ac.il (Dmitry Batenkov) 

^Current address: Department of Mathematics, Massachusetts Institute of Technology, 
Cambridge, MA 02139, USA. email: batenkov@mit.edu 

^The research leading to these results has received funding from the European Research 
Council under European Union’s Seventh Framework Program, ERC Grant agreement no. 
320649. 


Preprint submitted to Elsevier 


October 24, 2016 






This so-called confluent Prony problem and its numerous variations appear in 
signal processing, frequency estimation, exponential fitting. Fade approxima¬ 
tion, sparse polynomial interpolation, spectral edge detection, inverse moment 

§ )blems and recently in theory of super-resolution (see 0 , Si El, [II [II [II 
Him and references therein). We comment on the specific assumptions 
made in the above formulation in Subsection 11.31 below. 


Several solution methods for Prony systems have been proposed in the literature, 
starting with the classical Prony’s method and including algorithms such as MU¬ 
SIC/ESPRIT, Approximate Prony Method, matrix pencils, Total Least Squares, 




Variable Projections (VARPRO) or ^l minimization 
and references therein). While the majority of these algorithms perform well on 
simple and well-separated nodes, they are somewhat poorly adapted to handle 
either multiple/clustered nodes (the root extraction/eigenvalue computation be¬ 
coming ill-conditioned), large values of N (the quadratic cost function is highly 
non-convex w.r.t to the unknowns aij) or non-Gaussian noise. Despite this, 
our recent studies iid suggest that these problems are only partially due 
to the inherent sensitivity of the problem (i.e. problem conditioning). Gen¬ 
erally speaking, introduction of confluent (high-order) nodes into the model 
leads, in some cases, to improved estimation of the parameters - as indicated 
by the reduced condition number of the problem. In particular, we argue that 
while for NS ^ 1, where 6 is the node separation (see Definition [5] below), and 
D = (l,l,...,l) the existing methods might be close to optimal, there is a gap 
between theory and practice in the “near-collision” situation NS <C 1 and high 
multiplicity, even if the noise is sufficiently small. 


Decimation is a particular regularization for near-colliding Prony systems, which 
was first proposed in 0| and further analyzed in 0| . The essential idea is that if 
NS <Cl, then after taking the decimated sequences {mp, Wp, m 2 p, ■ ■ ■, ™(fi-i)p}, 
where p S N is not too large, as measurements, and solving the resulting square 
system, we can get accuracy improvement of the order of p~‘^^ for the node 
Zj - compared to the error in the case p = 1. Numerical studies carried out 
in 01 (see Subsection 12.21 below for further details) indicate that in this case, 
the best possible resulting accuracy is very close to the accuracy obtained by 
least squares - as quantified by the “near-collision condition number”. The work 
0 does not suggest any practical solution method which achieves the above 
accuracy, and in the present paper we propose to fill this gap. 


1.2. Our contribution 

In this paper we focus on developing an accurate solution method for Problem 
[T]in the near-collision regime NS 1, in the case of a single cluster. 

We propose a novel symbolic-numeric technique, “decimated homotopy”, for this 
task. The approach is an extension of the method used in 0 for the case s = 1 
(i.e.only one node), and its main ingredients are: 

1. decimating the measurements; 
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2. constructing a square polynomial system for the unknowns {zj} (in Q 
this was a single polynomial equation) ; 

3. solving the resulting well-conditioned system with high accuracy; 

4. pruning the spurious solutions and recovering the solution to the original 
system. 

Step 2 above is a purely symbolic computation based on the structure of the 
equations o, while for step 3 we chose the homotopy continuation method for 
polynomial systems, due to the fact that it will provably find the solution. We 
propose several alternatives for the pruning step 4, and discuss their efficiency. 
We also show that the proposed algorithm recovers the nodes with high accuracy 
- in fact, with near-optimal accuracy. Numerical simulations demonstrate that 
the algorithm is accurate as predicted, and outperforms ESPRIT in this setting. 
Some of the presented results have been published in abstract form in the pro¬ 
ceedings of the SNC’14 meeting Q, which took place in Shanghai, China, during 
July 2014. 

1.3. More on the assumptions 

Let us briefly comment on the specific assumptions made throughout the paper, 
and point out some current limitations. 

1. The setting \zj\ = 1 is common in applied harmonic analysis, where the 

prototype model for Problem [T] is to recover a Dirac measure / (x) = 
J2i=i ~ ^i) from the Fourier coefficients Ck (/) = ^ / (t) dt. 

Dropping this assumption will in general have severe consequences in terms 
of numerical stability of the problem. 

2. The confluent/high-multiplicity models are also quite common in inverse 
problems involving some sort of derivatives. In the most general formula¬ 
tion, the multiplicity structure D may be unknown, and it is an important 
question to determine it reliably from the data and other a-priori infor¬ 
mation. This is an ongoing research effort, and we hope that the methods 
of the present paper may serve as a building block for this goal. See also 
a discussion in Section |6] below. 

3. Similarly, extending the treatment to multiple clusters is an ongoing work. 
In this regard, our method can be regarded as a zooming technique. 

4. The noise level e is assumed to be sufficiently small. Quantification of 
the allowed noise level for the problems to be solvable is a closely related 
question and treated elsewhere, see e.g. [H,[^. 

1 . 4 . Organization of the paper 

In Section [5] we discuss in detail the relevant prior work, in particular accuracy 
bounds on (decimated) Prony systems [^, and the algebraic reconstruction 
method for the case s = 1 from The decimated homotopy algorithm is 

subsequently developed in Section |31 Analysis of the algorithm and its accuracy 
is presented in Section |4l Results of numerical experiments are described in 
Section [SJ while several future research directions are outlined in Section [SI 
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2. Accuracy of solving Prony systems, decimation and algebraic re¬ 
construction 

We start with brief recap of the Prony’s method in Subsection 12.11 Following 
in Subsection 12.21 we present numerical stability bounds, including in the 
decimated scenario, for the system ([T]). In Subsection 12.31 we discuss the “alge¬ 
braic” reconstruction algorithm for the system (H]) with s = 1, used in d t3|, 
and highlight some of its key properties, in particular the effect of decimation 
on its accuracy. 


2.1. Prony’s method 

The high degree of symmetry in the system of equations © allows to separate 
the problem into a linear and a nonlinear part. The basic observation (due to 
R. de Prony [^) is that the sequence of exact measurements {mk\ satisfies a 
linear recurrence relation 

d 

ruk+ici = 0, /c e N, (2) 

fco 

where {ci} are the coefficients of the Prony polynomial defined as 


Pix) := W{x-ZjY^ 

1=1 


d 


y^^cjx^. 


( 3 ) 


Thus, the system o can be solved for N = 2d hy the following steps. 


1 . 


Using (m, recover the coefficients {c^} of P (x) from a non-trivial vector 
in the nullspace of the Hankel matrix 

/mo • • • rhd-i rhd \ 


Hd-.= 

\md-i rhd 


m2d-i 


( 4 ) 


2. Recover the nodes {zj} by finding the roots of P (x) with appropriate 
multiplicities. 

3. Given the the nodes {zj}, recover the coefficients {aej} by solving a Van¬ 
dermonde linear system. 


2.2. Stability bounds and decimation 

Let us first introduce some notation. The number of unknown parameters is 
denoted hy R := d + s. 

Definition 1. The data space associated to the Prony problem is the a-priori 
set of possible solutions 

^ ■ — *1^ (^0,1) • ■ • ) Udi —1,15 ^1') • ■ • HQjS, ■ • ■ , Hds — l,s) € G . \Zj\ — ^ ^ 
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We also extensively use the notion of node separation, defined as follows. 

Definition 2. Let a; G S) be a data point as in ©• For i j, let Sij := 
\argZi — argZjl with the convention that Sij ^ tt. For i = 1,..., s the i-th node 
separation of x is 

niinJij. (6) 

In addition, we denote the global separation as 

5 = 5 (x) ■= min . 

i 


For any N ^ R, let the forward mapping Vj^f : S) —> C'^ be given by the 
measurements, i.e. for any a; G J) (see ([S])) we have 

V^r {x) ■■= ( too , • ■ , 


where to^ are given by (P). 


A standard measure of sensitivity 
terns is the following. 


iH 34 1 for well-conditioned polynomial sys- 


Definition 3. Let a; G 2) be a point in the data space. Assume that Jf^f {x) := 
dV^f (a;), the Jacobian matrix of the mapping V^f at the point x, has full rank. 
For a = l,2,...,i?, the component-wise condition number of parameter a at 
the data point x is the quantity 


N 

CN^,^{x) ■.= Y.\Jm^ (a^L,, 
2=1 


( 7 ) 


where is the Moore-Penrose pseudo-inverse of Jjg-. 

In @ we show that for N5 ^ 1, the Prony system ([T]) is well-conditioned 
as follows (up to changes of notation and a slightly different noise model, see 
footnote on page 4 in [6|): 

Theorem 1 (Theorem 2.1 in [il). Let a; G 2) he a data point, such that 5 = 
5 {x) > 0 and Od^-ij 0 for j = 1,..., s. Then 

1. The Jacobian matrix fjjg- {x) = dVj^ {x) G has full rank. 

2. There exist constants K, not depending on N and 5, such that for 

N>K-5-^: 


CN^.^n {x) 


< 


c(i). 


1 1 


( 8 ) 


It is easy to show (see e.g. Appendix A]) that the upper bound on CN^.^n 
is asymptotically tight. 
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On the other hand, as numerical experiments in Q show, when NS ^ 0 then 
the growth of is much more rapid than (obviously —>■ oo 

as the system becomes singular). As we now argue, this “phase transition” near 
NS ^ O (1) can be partially quantified by considering a sequence of decimated 
square systems. 

Fixing N = R, we have the following upper bound, which is tight. 

Theorem 2 (Theorem 2.2 in [^). Assume the conditions of Theorem [H and 
furthermore that N = R. Then there exists a constant not depending on 

X (and in particular on 5), such that: 

A natural question is whether increasing N can essentially improve the bound 
((9|) above. One possible answer is given by what we call “decimation”, as follows. 

Definition 4. Let p S N be a positive integer. The decimated Prony system 
with parameter p is given by 

s dj — l 

Uk := mpk = (aijp^) k^, fc = 0,1 ,..., i? - 1. (10) 

j=i e=o 


Definition 5. The decimated forward map —>■ is given by 

T’(P) (x) := (no ,.. , 

where a: £ J) is as in ® and Uk are given by (fTUl) . 


Definition 6. The decimated condition numbers CNff are defined as 

R 


(x) := 


i=l 


(*)} 


-1 


( 11 ) 


where (x) is the Jacobian of the decimated map (the definition applies 
at every point x where the Jacobian is non-degenerate). 


The usefulness of decimation becomes clear given the following result. 

Theorem 3 (Corollary 3.1 in @). Assume the conditions of Theorem [TJ As¬ 
sume further that NS* <ttR where S* := Sij (i.e. all nodes form a clus¬ 

ter). Then the condition numbers of the decimated system (1101) with parameter 
P* ■= Lf J satisfy 

i ( 12 ) 
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The intuition behind this result is that decimation with parameter p is in 
fact equivalent to applying the Prony mapping V-ji to a rescaled data point 
y:=TZp{x), where 


• ■ • ; —l.lj '2^1 5 • ■ • : ; ^ds — l,Si -^s) j ■- 

(6o,i, • ■ •, wi, ■ ■ ■, &o,s, ■ ■ ■, bd,-i,s, Ws)’^ = (13) 

(ao,i Odi-i,! • zf,..., ao,s • p°, ■ • ■, ad,-i,s ■ . 

Since for small 5^^'^ we have that mini^j \zf — Zj \ ~ (fl^ follows from the 

above and 

Experimental evidence suggests that decimation is nearly optimal in the “near¬ 
collision” region, i.e. 

CiVif) (a;) CN^^,n (x) , NS* < ttR. (14) 

We believe that it is an important question to provide a good quantification of 

m- 

From the practical perspective, the above results suggest a nearly-optimal (in 
the sense of conditioning) approach to numerically solving the system ([T]) when 
all nodes are clustered - namely, to pick up the R evenly spaced measurements 

|toO) Wp) • ■ •) i)p} 

and solve the resulting square system. 

An important feature of the decimation approach is that it introduces aliasing 
for the nodes - indeed, the system has Wj = Zj as the solution instead 
of Zj, and therefore after solving the algorithm must select the correct 

value for the root {wj)p . Thus, either the algorithm should start with an 
approximation of the correct value (and thus decimation will be used as a fine- 
tuning technique), or it should choose one among the p candidates via some 
pruning technique - for instance, by calculating the discrepancy with the other 
measurements, which were not originally utilized in the decimated calculation. 

2.3. Algebraic reconstruction 

Although many solution methods for the system (HD exist, as we mentioned they 
are not well-suited for dealing with multiple roots/eigenvalues. While averaging 
might work well in practice, it is difficult to analyze rigorously, and in particular 
to prove the resulting method’s rate of convergence. 

In we developed a method based on accurate solution of Prony system 

for resolving the Gibbs phenomenon, i.e. for accurate recovery of a piecewise- 
smooth function from its first N Fourier coefficients. This problem arises in 
spectral methods for numerical solutions of nonlinear PDFs with shock dis¬ 
continuities, and was first investigated by K.Eckhoff in the 90’s 0. The key 
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Algorithm 1 Algebraic reconstruction of a single node 

1. Set decimation parameter to p* := ■ 

2. Construct the polynomial Qp* (u) from the given perturbed measurements 

Qp, (u) := ^ (-1/ 

3. Set p to be the root of qp* closest to the unit circle in C. 

4. Choose the solution z* to m among the p* possible values of (p*)^ 
according to available a-priori approximation. 


problem was to develop a method which, given the left-hand side {mk}k=o 
m with error decaying as |Amfc| ^ fc ^ , would recover the nodes {zj} with 
accuracy not worse than \Azj \ ^ 

Our solution was based on two main ideas: 

1 . Due to the specifics of the problem, it was sufficient to provide a solution 
method as described above in the case of a single node, i.e. s = 1. 

2. The resulting system was solved by decimation, elimination of the linear 
variables {oij}, and polynomial root finding. 

The elimination step is a direct application of the recurrence relation ([2|) for the 
coefficients of the Prony polynomial (O, as follows. The (unperturbed) system 
o for s = 1 reads 

d-l 

ruk = z^''^aik^. ( 15 ) 

e=o 

The corresponding decimated system m with parameter p is 

d-l 

nk = rupk = {z’Pf ^ k = 1 ,2,.. . ,d + 1. 

1=0 

Denote p := Then clearly the sequence {uk} satisfies ''T‘k+eCi=0, where 

the Prony polynomial is just (x — p)'^ = That is, Cf=(—1)^ 

and we obtain that p is one of the roots of the unperturbed polynomial 

qp {u) := (- 1 )^ ( 16 ) 

e=o ^ 

The algebraic reconstruction method for s = 1 ([^, Algorithm 2]) is summarized 
in Algorithm [1] on this page. 

The key result of is that as A —>■ oo (and therefore p* —oo as well), and 
assuming perturbation of size e for the coefficients {m^}, all the d roots of qp* 
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remain simple and well-separated, while the corresponding perturbation of the 
root p* is bounded by 

\P-P\^ ^ |5 - z| < N-'^e. 


Thus, the method is optimal - recall the condition estimate ([S]). 

The pruning step 3 was shown to be valid since the unperturbed polynomial 
qp* has only one root on the unit circle. Regarding step 4, it was shown that 
a sufhcientR accurate initial approximation can be obtained by the previous 
method of Q. 

Remark 1. Decimation acts as a kind of regularization for the otherwise ill- 
conditioned multiple root. To see why, consider the case d = 2. Then we have 

TTik = (ao -I- kai). 


The Prony polynomial is P (x) = {x — zf ^ and thus for each fc G N the point z 
is a root of 

qf (u) := mkU^ - 2umk+i + mfc+ 2 - 
As fc —> c», the above polynomial “approaches in the limit” 



k 


—>■ {u 


Thus, a “non-decimated” analogue of Algorithm [T] (such as i[I3) would be 
recovering an “almost double” root z *, and it is well-known that the accuracy of 
reconstruction in this case is only of the order ^/s when the data is perturbed 
by e. On the other hand, qp (u) = mpU^ — 2um2p + rnsp, and as p —> c» it is 
easy to see that 


—>■ aip {u^ — Apu + 3p^) 


aip {u - p){u- 3p ), 


i.e. the limiting roots are well-conditioned. 


3. Decimated homotopy algorithm 

In this section we develop the decimated homotopy algorithm, which is a gen¬ 
eralization of Algorithm [T] to the case s > 1. We assume that the multiplicity 
D is known, and the noise level is small enough so that accurate recovery of the 
nodes by solving the decimated system (TTOll according to Theorem [3] is possible. 
Recall that the feasible solutions are restricted to the complex torus 

T^:={zGe: |(z)J = l, z = l,...,s}. 
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3.1. Construction of the system 

Consider the decimated system with fixed parameter p. Denote Wj := zj. 
The decimated measurements {nk}^^Q satisfy for each A: = 0,..., s — 1 

d 

^ ^ — 0 , 

i^O 

where Ci are the coefficients of the Prony polynomial 

s d 

-P (a^) = n 

3=1 fco 


Let fjiixi,... ,Xd) denote the elementary symmetric polynomial of order i in d 
variables. Then we have 


/ 1 \a—t. 

Cl = (— 1 ) (Jd-l 


Wi, 


,Wi, 


■.= Tl(wi,...,'Ws). (17) 


xdi 


y-ds 


Thus the point w = {wi ,..., rcg) G T® is a zero of the s x s polynomial system 

(18) 


|/fc^^ (“) := (ti) = o| 


k=0....,s-l 


This Hankel-type system is therefore our proposed generalization to the poly¬ 
nomial equation (1161) . 


Example 1. s = 2, dj 


fo (u) 
fl (u) 


riQvfu^ + ni 
niufu^ + 712 


= 2. The system (IT8)) reads 

—2u^U2 — 2niti|) + n2 (uf + 4nin2 + 'f*!) + (~2ni — 2u2) + 

—2ufu2 — 2niti|) + n3 (lij -|- 4tiin2 + 'f*!) + (—2ni — 2u2) + ns 


3.2. Recovering the solution 

Generalizing the root finding step of Algorithm [U we propose to use the ho- 
motopy continuation method in order to find all the isolated solutions of the 
(perturbed) system (IT^ . 

While a-priori it is not clear whether the variety defined by (fT51) has positive¬ 
dimensional components, we show in Section |T] below that our wanted solution 
is indeed isolated, and therefore the homotopy will find it. Furthermore, by 
analyzing the Jacobian of the polynomial map at the solution and comparing it 
with the estimate m, we show that the obtained accuracy is optimal. 

We now consider the question of how to recover the correct solution of the 
original problem ([T]) from among all the isolated solutions 

S = {ui,.. .,us} 

of (ITS)) . Two issues need to be addressed. 
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1. Spurious solutions. Transformation to the Hankel-type polynomial system 
introduces spurious solutions m G C® which are not in T® and therefore 
cannot be equal to the p-th power of a solution to the original system (HD- 

2. Aliasing. Given a solution u G 5 fl T® to the p-decimated system ([T51) . 
there are in general p® possible corresponding solutions to o, namely 

Zp (ti) := {(^1,..., G T® : (zif = (it) J . 

Overall, the set of all possible candidate solutions is 

^ U ■ 

uGS 

We suggest several pruning strategies. 


Exhaustive search. Since we have access to the original system, we can select 
the solution giving the smallest residual for the non-decimated equations (IT51) . 
i.e. 


d-l 

z* := argmin,^g5 ^ (z) 

fe =0 


(19) 


Optionally, the upper index in the summation in (1191) can be increased to TV — 
d- 1. 


Pre-filtering. Instead of considering the whole set Q, one can first prune S and 
work with the solution which has smallest distance to T®, i.e.: 


M* ^ argmin niax |l-|(Mfc)J 

2 = 1,...,S 






g^Zp (u*). 


( 20 ) 


Using an initial approximation. In some applications it may be possible to ob¬ 
tain an a-priori approximation to the location of the desired solution z. So 
suppose that the algorithm is provided with Zina G T® and a threshold p > 0, 
then additional pruning can be achieved by setting 

G^{z&g-. \z - Zinit\ ^ rj} . (21) 

The strategies can be mixed, i.e. pre-filtering (l20l) can be followed by either 
m or m etc. 

The decimated homotopy is summarized in Algorithm [5] on the next page. 

In the next section we prove that for small enough noise, the exhaustive search 
(1191) is guaranteed to produce the approximation to the original solution, which 
is near-optimal. However, this strategy quickly becomes prohibitive. The prun¬ 
ing strategy (|20l) combined with (1211) is empirically shown to be as accurate and 
much faster, see Section [5l 
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Algorithm 2 Decimated homotopy algorithm 
Given: (mo, ..., fhN-i) G C^, multiplicity D. 

1. Set decimation parameter p* = • 

2. Construct the system 

■^p* : ifk’^ (“) — ™p*(fe+i)'^i (“) = 0 

I i=0 

3. Solve T-Lp* by homotopy continuation method. Let S be the collection of 
its isolated solutions. 

4. Select z* &Q using either (HU) , (EnD, (EID, or a combination thereof. 



4. Analysis 

In this section we analyze the proposed algorithm in the case of a single clus¬ 
ter. We show that for small enough noise level, the exhaustive search (HU) is 
guaranteed to produce a near-optimal approximation to the original solution. 


^.1. Statement of the results 

Let a; G 2) be a data point ([^ satisfying the conditions of Theorem [31 and 
p* be the corresponding decimation parameter. The corresponding decimated 
measurements Uk = Uk (x) are defined as in m- 

s 

rik = ^ bijk^, 

3=1 t=D 


where Wj = Zj and bgj = a^jp*^. As before, w = w {x) = {wi,... ,Ws). 

Definition 7. Let Vx (t) G denote the Jacobian matrix of the system 

(fTSl) at the point u = t: 


Vx (<) := 


/ dfk (a;) 

V duj 


k=0,...,s-l 


Theorem 4. Let the conditions of Theorem\^be satisfied. Then detDx (w) ^ 0. 
Therefore, u = w is an isolated solution of m- 

We adopt the definitions from (in particular, see Section 3.2.3 and Section 
9.1.2) for measuring linearized sensitivity of the solutions of empirical polyno¬ 
mial systems with respect to perturbations of the coefficients. Not surprisingly, 
they parallel our earlier definitions of conditioning (see Section E]). 

Let j be a multi-index, and denote by the monomial . For k = 

0 ,..., s — 1, let ak 3 denote the coefficient of in the equation number k in 
m- Finally, let Jk denote the set of multi-indices j for which akj 0. 
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Definition 8. Let n be a zero of the system m- Assume that the coeffi¬ 
cient Ofc j is perturbed by at most Aakj- The linearized sensitivity of the i-th 
component of u to the change in the coefficients of the system is 


s-l 


Ki{u) I ^ |m^| I maxIAafcjl, 

\3^Jk 


k=0 


where Ki^k are the entries of the matrix {Dx (m)} ^. 

Theorem 5. Let the conditions of Theorem be satisfied. For small enough 
noise e 1 m the right-hand side of ([T]), Algorithm\^ with (1191) recovers the 
solution u — w of (I18|) with accuracy 

R-di 


Ki (w) ^ 


1 


1 


Consequently, the original solution z of o is recovered with accuracy 




R-di 


€. 


( 22 ) 


(23) 


Here and are constants depending only on the multiplicity vector D. 


Comparing the bounds (1231) and (TT^ , we conclude that the proposed algorithm 
is optimal, up to a constant, in the considered setting. 


4 . 2 . Proofs 

We start by deriving explicit expressions for the entries of V. 
Lemma 1. For j = 1,..., s and arbitrary fc £ N we have 


dfk 

dui 




Proof. Considering the coefficients as functions of re, we have the identity 

d 

fk ({rife (re)} , re) = ^ nk+^ (re) r, (re) = 0. 

2=0 


Thus, for each Wj the total derivative {{nk} , u) vanishes on u (re) = re. By 
the chain rule 


— fk({nk},u) = 
dwi 


dfk dfk du( 

^ duk-i-i dwj ^ dui dwj 

d 1 

, fc + 2—1 


= (w) (k + i) ^ ^ bi^j (k + 0^ + 1^ i'w) 


i=0 

0 . 


r=o 
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Let Tj (fc) denote the following polynomial in k of degree dj: 


Tj (fc) := 


dj-1 

e=o 


Then we have 


(w) = -w] ^ Y (^) (k + i). 


2=0 


( 24 ) 


We now employ standard tools from finite difference calculus [I3- Consider the 
right-hand side of (1241) as a discrete sequence depending on a running index k. 
Let E = Efc denote the discrete shift operator in k, i.e. for any discrete sequence 
g (fc) we have 

Eg{k) = (E5) (k) =g{k + l). 


Let us further denote by A := E —I the discrete differentiation operator (I is 
the identity operator). Now consider the difference operator 


£j (25) 
2=1 


Recall the definition of Ti from dni- Opening parenthesis, we obtain that for 
any g(k) 

d 

^j9 (^) = X! g{k-\-i). 

2=0 


Therefore 


dfk 

duj 


(w) 


-w^ (fc). 


Since the linear factors in (ESI) commute, we proceed as follows: 


-w) (k) = -w’^ ^Y[ {wjE-Wir)‘^' (k). (26) 

It is an easy fact (e.g. [13) that for any polynomial p (k) of degree n and leading 
coefficient ao, we have that 


A"p (fc) = aon!. 

Since rj (k) has degree dj, we obtain that 

A"V,(fc) = d,!5d,_i.,. (27) 

Eurthermore, applying the operator WjFi —Wil to a constant sequence c(fc)=c 
gives 

(wj E —Wi I) c = {wj — Wi) c. (28) 
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Plugging (l?fl) and (E51) into ((^ we get: 


dfk 

dui 


(w) = n 

i¥=j 


completing the proof of Lemma [T] 
Example 2. For s = 3, dj = 2 we have 


V = 


—(uii — W2)‘^ {wi — uis)^ —2bi2W2 (wi — 102 )^ («'2 — 

— 2 f)iiuij (wi — W2)^ (wi — ws)^ —2bi2'W2 (wi — W2)^ (j«2 — 'ii’3)^ 

2 f)iiuij (wi — W2)^ (wi — ws)^ —2bi2W2 (wi — W2)^ (w2 — 103)^ 


□ 


-2613103 (wi - 103)^ (W2 - 103)^ 
—2613103 (toi — 103)^ (102 — 103)^ 
—2613103 (101 — 103)^ (102 — 103)^ 


Definition 9. Let V (w) denote the s x s Vandermonde matrix on the nodes 
{wi,... ,Ws} ■ For example, if s = 3 we have: 


V(w) 


1 

1 

1 

Wi 

W2 

W3 

2 

2 

2 

Wi 

wi 

wi 


Corollary 1. Let y = TZp* (x) where TZp* is the scaling mapping (H. Denote 
by B (y) the following s x s diagonal matrix: 


B{y) := diag^.^i^ ,, 


dj ! -1 ,j {wj Wi ) 


Then we have the factorization 


(to) = V (w)B(y). 


(29) 


Proof Directly from Lemma [TJ 


□ 


Proof of Theorem^ According to our assumptions, both V (w) and B (y) are 
non-singular (in particular, since Udj-ij ^ 0 and mini^j \wi — Wj \ > 0). Using 
(l2^ we conclude that D is invertible. The conclusion that u = w is isolated is 
a standard fact about multivariate nonlinear systems, see e.g. 0. □ 


Proof of Theorem O From (|29|) we have 




diagi=i,...,s 


1 

dilbdi-ij 


e=/=i 



V(w)~^. 


15 







Let V (w) ^ s- Using the classical estimates by Gautschi [u 

we have that 


^ ki.fei < n i,„, i 


-1 


k=0 




Therefore 


\Wj -Wi\ 




i#* 


fc=0 


— 1 , 






2 S-1 


k=0 


di\ \adi-i,i\P 


*d 


— n “ Wi) 


-dt-l 


l^i 


di\ \adi-i,i\p*‘^'~^ \p*5^^ 


R—di — 1 


1 f 1 


R—di — 1 


dil \adi-i,i\P*^ ^ 

Now clearly there exists a constant Ci = Ci (s) for which 


( y^ [ic^ l maxIAafej l < Cie. 
jeJk ) ^ 


Since < 7 ri? and p* = 

Ki (w) ^ 


we have the bound 

C 2 - 

1 

( 1 



C 3 

1 

( 1 




R-di-1 

R-di 


€ 


proving with = C 3 . 

Because extraction of p*-th root reduces error by a factor of p*, this immediately 
implies dlSl) with C^) = i?cU). 

Since the homotopy algorithm converges to the exact solution of the approximate 
system dUD, and since e is assumed to be sufficiently small, the exhaustive search 
(fTOl) must produce the exact solution to the perturbed Problem [T]- otherwise 
there would be two non-proportional vectors in the nullspace of the Hankel 
matrix Hd and this is impossible since the rank of Hd is known to be 
exactly d. □ 


5. Numerical experiments 
5.1. Setup 

We chose the model ® with two closely spaced nodes, varying multiplicity and 
random linear coefficients ■ 


16 










We have implemented two pruning variants: 

1. the exhaustive search m-, 

2. combination of (OiH) and (I^Tj) (referred to as filtering below). 

Choosing the overall number of measurements to be relatively high (1000-4000), 
we varied the decimation parameter p and compared the reconstruction error for 
Algorithmic with filterin g ab ove (referred to as DH below), and the generalized 
ESPRIT algorithm (3.[^.135| (see also [1|), one of the best performing subspace 
methods for estimating parameters of the Prony systems o with white Gaus¬ 
sian noise Ck- The noise level in our experiments was relatively small. 

In addition to the reconstruction error, for each run we also computed both the 
full and decimated condition numbers and CN^f^ from their respective 

definitions ([7]) and CD- 
Additional implementation details: 

1. PHCPACK Release 2.3.96 was used as the homotopy continuation 
solver. It was called via its MATLAB interface PHCLab [^. All tests 
were run on Apple MacBook Pro 2.4 GHz Intel Gore i5 with 8GB RAM 
under OSX 10.10.4. 

2. We used the value p k, N~^ for the heuristic (ED- 

3. The node selection in generalized ESPRIT was done via fc-means clustering 
on the output of the eigenvalue step. 

5.2. Results 

The results of experiments are presented in Table [I] Table El and Figure [T] on 
pageim They can be summarized as follows: 

1. The exhaustive search is accurate, but time-prohibitive even for moderate 
values of p (the number of solutions considered is on the order of s!fi®p®). 

2. The accuracy of DH surpasses ESPRIT by several significant digits in the 
near-collision region N6 <C 1. 

3. DH achieves desired accuracy in larger number of cases. 

Some additional remarks: 

1. The number of solutions of the system (fTSl) was equal to s!d® (s=number 
of nodes, d = dj=degree). 

2. Running times are better for DH when p is large, because the selection 
step of Algorithm E] is 0{N), while the cost of full SVD for ESPRIT is 
O [N'^R) . See in particular Table El We did not observe any unusual 
memory consumption during the execution of the algorithms. 

3. Gondition number estimates are somewhat pessimistic, nevertheless in¬ 
dicating the order of error decay in a relatively accurate fashion. The 
periodic pattern is well-predicted by the theory, see Q. 
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Nodes=[1.3000 1.2950] 
Order=[2 2 ] 
MaxN=2000,SNR=157.0 


Nodes=[1.3000 1.2800] 
Order=[2 2 ] 
MaxN=2000,SNR=136.9 



Figure 1: DH vs. ESPRIT. Also plotted are condition numbers, both full and 
decimated, as well as the threshold ry used in the experiments. 











p 

time (sec) 

rec.error 

3 

0.9 

0.0002 

4 

1.1 

0.0001 

6 

1.1 

0.0004 

8 

1.0 

0.00008 


P 

time (sec) 

rec.error 

1 

1 

0.005 

3 

8 

0.001 

5 

25 

0.0002 

10 

99 

0.00009 


(a) Exhaustive search. Pruning step (b) Filtering. Overall timings 
only. 


Table 1: Running times for exhaustive search and filtering. The experimental 
setup is the same as in Figure fTbl Each table describes a separate experiment, 
therefore the reconstruction errors are slightly different. 


P 

ESPRIT 

DH Create 

DH Run 

DH Select 

120 

0.13 

0.84 

0.11 

0.001 

169 

0.21 

0.75 

0.09 

0.001 

238 

0.21 

0.74 

0.10 

0.002 

335 

0.34 

0.89 

0.10 

0.003 


Table 2: DH vs ESPRIT timings (sec). The columns for DH correspond to the 
system construction, PHC run and the pruning steps. 


6. Discussion 

In this paper we presented a novel algorithm, Decimated Homotopy, for numer¬ 
ical solution of systems of Prony type ([T]) with nodes on the unit circle, \zj \ = 1 
which are closely spaced. Analysis shows that the produced solutions have near- 
optimal numerical accuracy. Numerical experiments demonstrate that the prun¬ 
ing heuristics are efficient in practice, and the algorithm provides reconstruction 
accuracy several orders of magnitude better than the standard ESPRIT algo¬ 
rithm. The pruning (EUl) will be justified in the case that the system TLp* has no 
spurious solutions on the torus T'*. This seems to hold in practice, and therefore 
a theoretical analysis of this question would be desirable. On the other hand, 
initial approximations of order 77 Ri N~^ can presumably be obtained W existing 
methods with lower-order multiplicity (similar to what was done in [7]). 

Another important question of interest is robust detection of these near-singular 
situations, i.e. correct identification of the collision pattern D. While the inte¬ 
ger d can be estimated via numerical rank computation of the Hankel matrix Hd 
(HI) (see e.g. [13 and also a randomized approach [23 ); the determination of the 
individual components of I? is a more delicate task, which requires an accurate 
estimation of the distance from the data point to the nearest “pejorative” man¬ 
ifold of larger multiplicity, and comparing it with the a-priori bound e on the 
error. We hope that the present (and future) symbolic-numeric techniques such 
as [I1,|33, combined with description of singularities of the Prony mapping Vj^ 
13, will eventually provide a satisfactory answer to this question. 
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As we discuss in Q, decimation is related to other similar ideas in numerical 
analysis and signal processing 2 ^ 2^ 13 ■ symbolic-numeric literature 
connected with sparse numerical polynomial interpolation (i.e. in the noisy 
setting), the possible ill-conditioning of the Hankel matrices Hd can be overcome 
either by random sampling of the nodes {zj} [S 
introduced affine sub-sequence approach |2 
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I or by the recently 
for outlier detection (see also 0), 


which is in many ways similar to decimation. 


It would be interesting to establish more precise connections of our method to 
these works. 
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